Numerical solution of Riemann-Hilbert problems: 
random matrix theory and orthogonal polynomials 

Sheehan Olver and Thomas Trogdon 
October 9, 2012 



Abstract 



(N 

O 

(N 

o 

o 

oo 



B 



> 

On 
On 

d 
X 



In recent developments, a general approach for solving Riemann-Hilbert problems 
numerically has been developed. We review this numerical framework, and apply it to 
the calculation of orthogonal polynomials on the real line. Combining this numerical 
algorithm with an approach to compute Fredholm determinants, we are able to calcu- 
late level densities and gap statistics for general finite-dimensional unitary ensembles. 
We also include a description of how to compute the Hastings-McLeod solution of the 
homogeneous Painleve II equation. 

1 Introduction 

We are concerned with calculating random matrix statistics for Hermitian invariant ensem- 
bles; i.e., n X n random matrices 
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whose entries are distributed according to 



^(n-l)n - i^(n-l)n 



Ml + iMl\ 
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where Z„ is the normalization constant and 



dM = J]dM,,J](dM5dM^.) 



1=1 



i<j 



The eigenvalue statistics of invariant ensembles are expressible in terms of the kernel 

^^_^e-^mvi.)+viy)) vr„(a;)7r„_i(y) - 7r„_i(x)7r„(y) 



27ri X — y 

where iTk are the orthonormal polynomials ttq, tti, . . . with respect to the weight 

e-"^(") dx, 
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and ^ 

7„_i = 27ri j ']Tn-i{x)w{x) dx 

is a normalization constant. Particular statistics include the level density 

d/i„ = /C„(x, x) dx, 
describing the global distribution of eigenvalues, and the gap statistic 

det(/ - /C„|l2p), 

where det denotes a Fredholm determinant, describing the local distribution of eigenvalues; 
namely, the probability that no eigenvalue is inside the set VL. 

Gap statistics for invariant ensembles follow two principles of universality. For x in the 
hulk — i.e., inside the support of the equilibrium measure — the gap statistic of a properly 
scaled neighbourhood of x will approach the sine kernel distribution: 

det(/-5i2(_ )) for S = . 

^ ' X — y 

This was proved rigorously in [6] by expressing the orthogonal polynomials in terms of a 
Riemann-Hilbert problem, so that asymptotics of 7r„ were determinable via nonlinear steepest 
descent. Moreover, the edge statistic — i.e., a properly scaled neighbourhood of oo — 
generically approaches the Tracy-Widom distribution: 

det(/-.4U,,_,) for ^, AiWAi-(i/)-Ai-MAifa) ^ 

X — y 

Underlying these two universality laws are Painleve transcendents; in the case of the Tracy- 
Widom distribution it is the Hastings-McLeod solution to Painleve II [16], whereas the sine 
kernel distribution is expressible in terms of a solution to Painleve V pTZ] • See Section |A] for 
a discussion of a numerical Riemann-Hilbert approach for computing the Hastings-McLeod 
solution of Painleve II. 

The statistics differ from universality laws for finite n and are no longer expressible 
in terms of Painleve transcendents. Hence our aim is to calculate the finite-dimensional 
statistics to explore the manner in which the onset of universality depends on the potential V. 
To accomplish this task, we will calculate the associated orthogonal polynomials numerically, 
also using their Riemann-Hilbert representation, via the framework of [2T1 [22] • By deforming 
the contours appropriately, we will achieve a numerical method that is uniformly accurate 
for large and small n, as shown in [23] . 

We will see in our numerical experiments that the onset of universality depends strongly 
on the magnitude of the equilibrium measure: where eigenvalue density is small, finite n 
statistics differ from universality behaviour greatly. 

We begin with a demonstration of the numerical calculated finite-dimensional random 
matrix statistics (Section |2|. Importantly, because we do not require the knowledge of local 
parametrices, our numerical approach continues to work for degenerate potentials, such as 
those that arise in the study of higher order Tracy-Widom distributions [3] . We describe the 
manner in which orthogonal polynomials can be reduced to a Riemann-Hilbert problem that 
is suitable for numerics (Section [s]). We then review the numerical method for Riemann- 
Hilbert problems (Section |4]), based on the deformations of [6j. This includes the result 
that the numerical approximation is uniformly accurate when the contours are appropriately 
deformed (Section [s]), without the use of classical local parametrices. 
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Figure 1: Calculated level densities for the GUE for n = 3, 10 and 100, compared to 
histograms. 

Remark An alternative to the approach advocated in this paper is to calculate the 
orthogonal polynomials directly for each n via Gram-Schmidt and numerical quadrature. 
For small n, this is likely to be more efficient. However, it is well known to be prone to 
instability [E]; moreover, the calculation must be restarted for each n as the weight e""^ 
changes. On the other hand, the RH approach has computational cost independent of n, 
making it more practical for investigating large n behaviour. 



2 Random matrix theory 

Recall that 



27ri X — y 



and 



}Cn{x,x) = — «(x)7r„_i(x) - <_i(a;)7r„(x)). 



In this section, we use the approach of numerically calculating vr^ and 7„_i7r„_i that we de- 
velop below to compute the finite n statistics. What will be apparent in the numerical results 
is that the behaviour of local statistics is tied strongly to the global density of eigenvalues 
near the region; i.e., the magnitude of the level density. 

For unitary invariant ensembles, the level density is the distribution of the counting 
measure. This is precisely 

]Cn{x,x) 

n 

In Figure [l| we compare the (numerically calculated) GUE (i.e., V{x) = x'^) level density 
for n = 3, 10 and 100 to a histogram, demonstrating the accuracy of the approximation. 
(Because the polynomials involved are Hermite polynomials, we can also verify the accuracy 
directly.) This shows the standard phenomena that the distribution exhibits n "bumps" of 
increased density, corresponding to the positions of the finite charge energy minimization 
equilibrium; i.e., the Fekete points. 

In Figure [2| we plot the finite n level densities for the potential 

T.. A 4 3 8 

V(x) = X H h -X, 

^ ' 5 15 20 5 ' 

which is an example of a potential whose equilibrium measure vanishes at an endpoint, and 
hence the edge statistics follow the higher order Tracy-Widom distribution Interestingly, 
this change in edge statistic behaviour is not just present in the local statistics, but clearly 
visible in the decay of the tail of the global statistics. 
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Figure 3: The calculated probability that there are no eigenvalues in (— s, s) for the GUE 
(plain) versus Monte Carlo simulation (dashed), for n = 50 (left) and n = 100 (right). 

We now turn our attention to local gap statistics, which are described by the Fredholm 
determinant 

det(/ - /C„|l2p). 

Using the method of Bornemann [3], we can calculate the determinant, provided that the 
kernel itself can be evaluated. Thus, we can successfully calculate finite gap statistics by 
using the RH approach to calculate 7r„ and 'jn-i'^n-i- In Figure [3} we plot the gap statistics 
versus a histogram for the GUE in the interval (— s, s). 

To see universality in the bulk, we have to scale the interval with n; in particular, we 
need to look at the gap probability for 



Alternatively, /C be replaced by its asymptotic distribution to get 

" = - + ^- 

mp[x) 

where d/i = iIj{x) dx is the equilibrium measure of V. For x inside the support of fi, this 
statistic approaches the sine kernel distribution. We demonstrate this in Figure |4] for the 
degenerate potential, showing that the rate in which the statistics approach universality 
depends on the magnitude of the equilibrium measure. 

We now turn our attention to edge statistics. In the generic position (i.e., when the 
equilibrium measure has precisely square root decay at its right endpoint), the gap probabihty 
for 
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Figure 4: The calculated probability that there are no eigenvalues in the scaled neighbour- 
hood X + ^^^l.-^ for n = 50, 100,200 and 250 for a; = 1 (left) and x = 1.5 (right), for the 

potential V{x) = ^ — j^x^ + §5 + 5^- 




Figure 5: The equilibrium measure for V{x) = ~ x (left) and the scaled gap statistic for 
n = 10,20,40 and 80 (right). The dashed line is the Tracy-Widom distribution (n = 00). 



tends to the Tracy-Widom distribution; here c is a constant associated with the equilibrium 



measure, see Section 3.2 for its precise definition and the numerical method for its calculation. 



In Figure M we plot the computed equilibrium measure for V{x) = — x (computed as 



described in Section 3.1), and its scaled edge statistic for increasing values of n. While 
the finite statistics are clearly converging to the Tracy-Widom distribution, the rate of 
convergence is much slower than the convergence of bulk statistics where the density of the 
equilibrium measure is large. 

Remark There are several methods for calculating universality laws — i.e., n = 00 statis- 
tics — including using their Painleve transcendent representations, see [2] for an overview. 
An additional approach based on RH problems is to represent, say, 

9slogdet(/ - S\l2{-^s,s)) 

as a RH problem. This can be solved numerically for multiple choices of s, and the results 
integrated numerically, see [5| for examples in the degenerate case. This will be accurate 
in the tails, whereas the Fredholm determinant representation that we use only achieves 
absolute accuracy. However, we are not aware of similar RH problems for finite n. 
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3 Orthogonal polynomials 

We wish to calculate the monic polynomials vro(x), 7ri(x), . . . with respect to the measure 

e-"^(") dx 

supported on the real line. Consider the following RH problem: 
Problem 1 [13] The function 

7r„(z) C[7r„e-"^](;2) 



where 



1 



7„_i = 2m 

solves the RH problem 



TTn-l{x)'w{x) dx 



r, = yi( , and Y ' 
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To apply the numerical method described in Section]!], we must transform the RH problem 
for Y into a suitable form for numerical solution. To accomplish this, we will transform Y 
by representing it explicitly in terms of new functions which satisfy the following properties: 

1. y I— 7- T so that T ~ / at infinity. 

2. T i-> S" so that the oscillatory jumps of T become exponential decaying jumps of S. 

3. 5 I— 7- $ so that the jumps of $ are localized and scaled. 



3.1 Equilibrium measures 

Our first task is to remove the growth in y at oo. To accomplish this, we must compute a 
so-called ^^-function associated with the equilibrium measure of V: 

Definition 1 The equilibrium measure /i is the minimizer of 

jj 1^ 1 ^1 '^/^(^) + j V{x) d/i(x). 

In this section, we assume that the equilibrium measure of V is supported on a single interval 
(a, 6); a sufficient condition is that V is convex ]6j. (We remark that the below procedure 
was adapted to the multiple interval case in [19], and adapting our numerical procedure for 
computing orthogonal polynomials, and thence invariant ensemble statistics, to such cases 
would be straightforward.) 

With the correct choice of there exists g analytic off (a, 6) satisfying 

g+{x) + g-{x) = V{x) — i for a < x < b and g{z) ~ logz. 

The derivative (p = g' is analytic off (a, b) and satisfies 

4>+{x) + (j^^ix) = V'{x) for a < X < b and (l){z) ~ -. 
Given a candidate (a, b), we can describe all satisfying this property 
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Theorem 1 11 Denote the affine map from {a,b) to (—1, 1) as 

2z — a — b 



b — a 



Suppose we have 



..[X) 

where Tk is the kth order Chebyshev polynomial of the second kind. If 

+ (f>-{x) = V'{x) for X G {a,b) and 0(oo) = 0, 
then there exists a x such that 

°° Vh , ^ , . b — a b — a 



^ ^ 2y z — b^z — a 2y z — b^f z — 



for the inverse Joukowski transform 



J^^{z) = z - Vz - Wz + 1. 
Sketch of Proof This theorem foUows from Plemelj's lemma and the fact that 

J7\x)'' + J7\x)-^ 



where 



J, ^(x) = X — iVl — x\/l + X = hm J,^(x + ie) 



To achieve the desired properties, we want (p to be bounded: 

Vo = and x = 0- 

We also want (hiz) ~ -: 



= 1- 



These two conditions give us a function 

F{a,b) : 



Vo 

{b - a)Vi - 8 



□ 



for which we want to find a root. We can calculate Vq and Vi to high accuracy using the 
trapezium rule applied to 

This calculation is trivially differentiable with respect to a and 6, hence we can easily apply 
Newton iteration to find a root of F. Convexity ensures that this root is unique . 

Once (a, b) are computed, we calculate 0(z) by using the discrete cosine transform to 
calculate the Chebyshev coefficients of V . We then have the equilibrium measure 

d/i = ^ - dx = ^ ~ ^^"''^^""^ VkUu-i{M^a,b){x)) dx 

fc=i 
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where Uk are the Chebyshev polynomials of the second kind. 
To calculate g, we compute an indefinite integral of |T9] : 



[z) dz 



b — a 
4 



^1 ( TT^ l0g^+ {M(a,b){^)) 



k=2 



fc-i 



k + 1 



k-1 



This formula was derived by mapping J^^{M(^a,b){z)) back to the unit circle, where it became 
a trivially integrable Laurent series. Note that g has a branch cut along (— oo, a) on which 
it satisfies: 

- g-ix) = 27ri. 
Choosing (arbitrarily) x G {a,b), we calculate 

i = V{x) - g+{x) - g-ix). 

The numerically calculated g consists of approximating Vk using the discrete Cosine 
transform and truncating the sum. Due to analyticity, the errors in these computed coeffi- 
cients are negligible, and the approximation of g is uniformly accurate in the complex plane. 
Hereafter, we treat the numerical g and the true g as equal. 



3.2 Scaling constant for edge statistics 

Associated with the equilibrium measure are the Mhaskar-Rakhmanov-Saff numbers. We re- 
express the constant as stated in [7] in terms of constants that we have already calculated: 
the support of the equilibrium measure and its Chebyshev coefficients. The equilibrium 
measure for the scaled potential V{M^J^^{x)) has support (—1, 1). Its equilibrium measure 
is 



M[,\)i^)m[a\)i^))dx 



b — a 



i^{M^a\){^))dx={b-a 



x^ 



An 



^ Vfcf/fc_i(x) dx 



k=l 



2n 



-h{x) dx, 



as in ^ (3.3)]. We define the constant 

■/^(l)2^^/=^ 



a 



k=l 



2/3 



as in |l7i (3.10)]. The scaling constant is thus 

[b-a) 



b — a 



-1/3 



.k=l 



2/3 



3.3 Lensing the RH problem 

We can now rewrite Y to normalize the behaviour at infinity: 



Y{z) 
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e 2 



„ 111 
e 2 



T(z) 
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Figure 6: The jumps of T. 



so that T I and has a branch cut along the unit interval, on which it satisfies 

' ( -| pn(g++g_+^-y)\ 

X < a OX X >h 
a < X < b. 



gi(9+-S-) 

1 

n{g--g+) ]^ 



We appeal to properties of equilibrium measures (see [26]) to assert that 

g+{x) + g-{x) + £ - y < 

for X < a and x > b, thus those contributions of the jump matrix are isolated around a 
and b. On the other hand, (?+ — g- is imaginary between a and 6, hence e^"^^+~^~'' becomes 
increasingly oscillatory on (a, b). We wish to deform the RH problem into the complex plane 
to convert oscillations into exponential decay. To accomplish this, we introduce the lensing 
as in Figure [61 where we rewrite T by 



T{z) = S{z) { 



^n{V^e-2g) I 
1 

^n{V-e^2g) I 



2 G S+ 
otherwise. 



By substituting 



9+ = V -g-- 
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we see that the oscillations have been removed completely from supp /x: 



Ma- -9+) 



M9+-9-) 



_1 Q'n{9+-9-) 

1 

-1 



MV-i-2g-) lM_l QniV-e-2g^) 



However, we have introduced new jumps on F-i- and T^, on which 

1 



S,=T,=T_ = 



MV-e-2g) I 



3.4 Removing the connecting jump 

We have successfully converted oscillations to exponential decay. However, to maintain 
accuracy of the numerical algorithm for large n, we must isolate the jumps to neighbourhoods 
of the endpoints a and b. Thus we require a function which satisfies the following RH problem: 



N+{x) = N^{x) 
The solution is [6] 

Niz] 



1 



for 



a < X < b and A^(oo) = /. 



1 iA ujz) 

-i I 2 



1 -i 
i 1 



2u{z 

i.e., is a solution to 

= iz/„(x) for a < X < b 



for uiz) 



z — a 



1/4 



and 



z/ 00 



1. 



An issue with using as a parametrix is that it introduces singularities at a and b, 
hence we also introduce local parametrices to avoid these singularities. In the event that the 
equilibrium measure 'il}{x) has exactly square root decay at the edges, asymptotically accurate 
local parametrices are known. However, if the equilibrium measure has higher order decay 
(d la the higher-order Tracy-Widom distributions the asymptotically accurate local 
parametrices are only known in terms of a RH problem. 

For numerical purposes, however, we do not need the parametrix to be asymptotically 
accurate: we achieve asymptotic accuracy by scaling the contours. Thus we introduce the 
trivially constructed local parametrices which satisfy the jumps of S in neighbourhoods of a 
and b: 



Pa{z) 




^ < axg{z — a) < 71 



-71 < argU — a) < 



I < arg(2; — a) < 



otherwise 



MV-i-2g) 



-n{V-£-2g) 
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Figure 7: The jumps of 



and 




2-K 



< arg{z ~ b) < TV 



-TT < arg(z — b) < 



2tt 



< arg{z - 6) < 



otherwise 



MV-i-2g) 



-n{V-l-2g) 



We can now write 



S{z) 




\z — a\ > r and \z — b\ > r 
\z — b\ < r 
\z — a\ < r 



The final RH problem for $ satisfies the jumps depicted in Figure [7j 

In practice, we do not use infinite contours. We truncate contours when the jump matrix 
is, to machine precision, the identity matrix. In all cases we consider here, after proper 
deformations the jump matrices are C°° smooth and are exponentially decaying to the iden- 
tity matrix for large z. The truncation of contours can be rigorously justified by solving a 
'nearby' RH problem with truncated contours. For a full discussion of this see Section [5] 
and, in particular. Lemma |4] below. We can then deform the remaining contours to be line 
segments connecting their endpoints. The resulting jump contour consists only of affine 
transformations of the unit interval. 

4 Numerical solution of Riemann— Hilbert problems 

We have reduced the orthogonal polynomial RH problem to the form 

$+(^) = ^~{z)G{z), zeT, $(00) = /, (1) 
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where F = U ■ ■ ■ U F"^ is a union of contours that are affine transformations of the unit 
intervah i.e., M*(F*) = (—1, 1) for an affine transformation M*. We use the notation [G; F] 
to refer to this RH problem and $ = [G; F] when $ is the unique solution. Assume this 
solution has the form 

$(2) = J + CrU, 

for some smooth function U. The RH problem ([T]) is converted into an equivalent singular 
integral equation (SIE) by substituting this assumed form into ([!]): 

I + C^U = {I + C^U)G. (2) 

We use the operator identity [B] 
to rewrite ^ 

U-C^U{G-I) = G-I. (3) 

It is well-known that the operators Cp are bounded from L'^{T) to itself for every contour 
we consider here. We use the notation 

C[G; r]U = U ~Cj^U{G - I), 

which is a well-defined bounded linear operator on L'^(T) provided G G L°^(F). 

Our numerical scheme consists of approximating f/ by a finite-dimensional sum of mapped 
Chebyshev polynomials. In other words, for x G F we approximate U{x) ~ Um{x), where 
we define 

Um{x) = U\n{x) = Uin{M\x)) for x G V\ 

fc=0 

for as-of-yet unknown coefficients f/^ G C^^^. If we are given the coefficients, we can evaluate 
C [G; F] Um pointwise by using an exact expression for the Cauchy transform of our basis: 

Proposition 1 j2Uf 

Cr.[TfcoM^](z) = C(_i,i)Tfc(M^(z)) 
Sketch of Proof Follows from Plemelj's lemma: 

C:(_i,i)Tfc(MXoo)) = C:(_i,i)Tfc(oo) = 

and 



'-(-1,1) '-(-1,1) 



TkiM\x)) = Tk{M\x)). 

□ 



Theorem 2 l2^\2Tj Define 



2 
ivr 



arctanh z 

^l+2L-fc/2J+fc / 1^ 1 

(l-22)(l+2[-fc/2j) 2^iy3^l-k/2p ^■ 



z2_l J 



2;^ (arctanh 2; — arctanh z ^) 

^fc-l-2L(fc+l)/2J / 

+ (i-^-2)(i+2L(fc+i)/2j) 2^iy 



1, 1 



for k = 
for k < 



+ L{fc+1)/2J ' 2-2-1 



for k > 0, 
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where 2-F1 is the hypergeometric function fJE/- Then 
where 

J-^{z) = z - Vz - iVz + 1 
is again an inverse of the Joukowsky transform J{z) = ^{z + z~^). 



Sketch of Proof This also follows from Plemelj's lemma after mapping to the unit circle: 

TkiJiz)) 



z^ + z ^ 



2 ' 

and relating the Taylor series of arctanh z to the hypergeometric function. 



□ 



We now choose the coefficients f/^ by enforcing (|3| to hold pointwise at a collection of 
= \m\ = m} + ■ ■ ■ + points. In other words, we choose points {2;^, . . . , z'^^_^ lying on 
each Fj, and solve the AN x AN linear system 

C[G-V]Um{z\,) = G{zl)-I. (4) 

We choose mapped Chebyshev points for the points: 

{z\,...,zl^} = !^M^-\-l),M^~' (^cosn 1-;^ M^-\l)| . 

This means every junction point of F is included in the collocation system, with multiplicity 
the number of contours emanating from the junction point. We denote these repeated points 
by 

{e + 0e^^S...,e + 0e'^-} 

where 61, . . . ,6l are the angles in which the components of F that include ^ as a junction 
point emanate from ^. But, as seen in Theorem [2| the Cauchy transform for our basis blows 
up at such points! To overcome this discrepancy, we assume that the solution satisfies the 
zero sum condition: 

Definition 2 Um satisfies the zero sum condition if, at every junction point (, it satisfies 

where Pi = —1 if the left endpoint ofV^ is Q, Pi = 1 if the right endpoint ofV^ is ( and Pi = 
if ( is not an endpoint of F* . 

We can define an alternate expression for the Cauchy transform at the junction points: 

Definition 3 For z not an endpoint ofT^, 

Cr4TkoM']{z)=Cr4TkoM']{z). 

Otherwise, for z^ the left endpoint of F* and z\^ the right endpoint, if 6 ^ 6i define 

Cr« [Tfc o M'] (4 + Oe'^) = + ir^ arg(-e'(^-^')) 
Cr.[n o M^](4 + Oe'^) = + ir,^ arg(e^(^-^')) 
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for 



a\ = (-1)'^ + ^ + /i.(-l)] + r]: log 



2m m 



27ri 



^ - ^°^%-[/i.-i(l)+/Xfc(l)]+r^og 



27ri iTT 

where 

z 



27ri' 



I ^ I 



When 9 = 9i define 

t/ie appropriate limits. 

The usefulness of this alternative definition is that it is equivalent to the standard Cauchy 
transform for functions which satisfy the zero sum condition: 



Lemma 1 If Um satisfies the zero sum condition, then 

CUm{z)=CUm{z). 

Sketch of Proof Let C be a junction point of F. From the asymptotic behaviour of arctanh z, 
we see near ( that if ( is an endpoint of F* , 

Cr4Tk o M^]{z) r^-^^og\z- 4 1 + CI, 

where 6 = aig{z — () and pi is defined as in Definition [2] If ( is not a junction point of F\ 
then 

Cr.[Tfc o M%z) ~ Cr^Tk o M^Q =■ CI,. 
(where C^^i. is 6* independent). Thus, 



i k 

= E E 



27ri 

i k i i k 



i k 

since the zero sum condition ensures that 

i 

The remaining constant CI which we refer to as the finite part, are precisely the constants 
we defined above. □ 

Thus, assuming the coefficients Ul are in the space so that Um satisfies the zero sum 
condition, we can replace ^ by 

C[C,T]Um{zi) = G{zi)-I. (5) 

This is justified by the following: 
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Lemma 2 If the linear system ^ is nonsingular, then the calculated Um satisfy the 
zero sum condition. 

Sketch of Proof Let C be a junction point, and assume for simplicity that = 1 or and 
Fj are ordered by increasing arguments 9i. Define 



$f =/ + C^f/m(C + Oe 



and define 



Gi = G{C + Oe'^0 

(i.e., the hmit of the jump along F*). The collocation system imposes 

But the definition of C imposes that 

= + (^^.+1 - 0^)S and $r = <f t + (^^1 + 27r - 

for 



2m 



These equations give 

=^IG,---Gl + S 



+ 271- 9l)G^ ■■■Gl + - d-^-l)G: ■■■Gl 



i=2 



+ 271- 9l)I + ^(0, - ---Gl 



i=2 



where we use the well-posedness of the RH problem: 

Gl ■ ■ ■ Gl = I- 

If 

L 

{9, + 271- 9l)I + Y,^9, - ■■■Gl 



i=2 



is nonsingular (the nonsingular junction condition) , then 5 = 0, implying the zero sum 
condition. 

Now suppose the nonsingular junction condition is not satisfied, and we replace the 
condition in the collocation system that 



with S = 0. We then have 



$,+1 = and $1 = 



Thus 



^^Gl =^t_,GL-lGL 



^L-l^l ' ' ' Gl-iGl 



Li 
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and the removed condition is still satisfied. In other words, the two linear systems are 
equivalent. 

□ 

In conclusion, the fact that the linear system is nonsingular implies that the numerically 
constructed 

^m{z) = I + CUm{z) 

is analytic off F and satisfies the correct jumps at the collocation points. We can thus 
recover approximations to orthogonal polynomials from ^rn by undoing the transformations 

r ^ T ^ 5 ^ $. 

We have one last task: we need to scale the contours so that the numerical algorithm 
remains accurate for all choices of n. 

4.1 Spaces 

We follow [23] and interpret the operator defined by applying Cr and sampling the resulting 
function at {zl} as mapping of piecewise polynomials to piecewise polynomials. The sampled 
function at the values {zl} can be identified with its unique piecewise-polynomial interpolant 
and we use Irn to denote this interpolation operator. Define L^^(r*) to be the space of 
matrices with entries being m^th order polynomials. When F = U ■ ■ ■ U has intersection 
points we define 

i=l 

We define L|^^(F) to be the closed subspace of L|^(F) consisting of functions that satisfy 
the zero sum condition. Thus XmCr is a well-defined linear operator from L|^^(F) to 
L^YiiV). As is mentioned in [23], XmC[G; F] maps to a proper subspace of L^(F). 

For each component contour F* and k G N"^ we define H^iJ'^) and Vr'^'°°(F*) in the usual 
way [23]. For the contour F 

L L 

i=l i=l 

Define H^(T) to be the close subspace H^(T) consisting of functions whose {k — l)th-order 
derivatives each satisfy the zero sum condition. Finally, for a Banach spaces X we use C{X) 
to denote the Banach space of operators on X with the induced operator norm. 

5 Uniform approximation 

In this section we describe how the convergence of the numerical approximation of RH prob- 
lems can be made uniform in a parameter. We also refer to this uniformity as asymptotic 
stability of the numerical method. In the case of orthogonal polynomials, the relevant pa- 
rameter is n, the degree of the polynomial. We refer to the results of [23]. We assume that 
we have a sequence of RH problems F„] depending on the parameter n. The theory of 
[23] requires some assumptions on F„. Assume 
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Figure 8: The pre-scaled used for non-degenerate endpoints and the pre-scaled fl^ used 
for first order degenerate endpoints. 

where {fi^}^-^^ are mutually disjoint and have the form 

ni = aiW + PI 

We assume r„ is the disjoint union of contours, each of which is an affine transformation of 
a fixed contour. Once we have this separation of r„ we attempt to solve the RH problem 
[Gn] r„] in an iterative way. We define the restricted jumps = Gnl^^ and the jumps after 
variable change Hl{k) = G{^{ai^k + (31), k E Qj. The following result can be found in [23] . 
For notational simplicity we suppress the dependence on n. But it is important to note in 
the general case every function, but no domain, will depend on n. 

Lemma 3 (Scaled and shifted solver) Assume $1 = [H^; Qi] and define $1 = $1 (^^^^r- j ■ 
Furthermore, for each j = 2, . . . ,1 define = ^i{a^k + and set 

^3 = [^3-1. ■ ■ ■ '^i.H.'^i', ■ ■ ■ '^7-1.' ' '^.•(^) = (^) • 

Then $ = $^ . . . solves [G„; r„] . 

This lemma states that we can treat each disjoint contour separately. We first solve a RH 
problem on one contour and modify the remaining jumps with the solution. This process is 
repeated until all contours are taken account of. 

We use the following rule of thumb to determine the proper scalings a^^: 

Assumption 1 // the jump matrix G has a factor e"^ and corresponds to a qth order 
stationary point (i.e., 6{z) ^ C{z — P^Y), then the scaling which achieves asymptotic stability 
is (a constant multiple of)al = n~^/'^ . 

In the case of a non-degenerate equilibrium measure, g{z) ~ Ca{z — a)^/^ and g{z) ~ 
Cb{z — 6)^/^; we thus scale like n~^/^: 

Qi = + a and ^1 = n-^'^^f + 6, 
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where n° is depicted in Figure [sj and the angle of the contours are chosen to match the 
direction of steepest descent. In the first order degenerate case (eg., V{x) = ^ — + to 
g{z) ~ Cb{z — hy/'^ and so we scale like n""^/^ at the degenerate endpoint: 

fii = n-^'^n^ + a and ^1 = n-'^/^Q^ + b, 

where Q,^ is depicted in Figure [s] (the angle is sharper to attach to the new direction of 
steepest descent). Higher order degenerate equilibrium measures will require higher order 
scalings, but this can be determined systematically by investigating the number of vanishing 
derivatives of the equilibrium measure. 

This is the final form of the RH problem that we used in the numerical calculations of 
Section [2j The remainder of the paper is concerned with proving that this scaled and shifted 
RH problem achieves asymptotic accuracy, i.e., the error does not grow as n becomes large. 

5.1 Conditions for uniform approximation 

A significant question is whether each of these smaller RH problems is solvable. From a 
practical numerical standpoint this possible issue does not seem to affect the conditioning of 
the method. From a theoretical standpoint this question is settled for large n in ^3J provided 
— >■ for all j as n ^ oo with some mild restrictions on 

Assumption 2 Assume that the jump matrix G is C°° when restricted to each component 
F* of F and decays to the identity matrix faster than any polynomial at each isolated endpoint 
of F and at oo oo G F. 

This is true in all cases we consider here. The following lemma is proved in [23] . 

Lemma 4 (Contour truncation) For every e > there exists an matrix-valued function 
and a hounded contour F^ such that 

• G^ = I onV \ T ^, 

• \\Ge — G'||i2(r)nL°°(r) < o,nd 

• \\C[G]V]- C[Ge]V^]\\c{L^(j)) < e||Cp ||£(L2(r)). 

Note that when the jump matrix G is the identity matrix then the solution of the RH 
problem is analytic across the jump. In practice we truncate infinite contours to finite 
contours when the jump matrix is within machine precision of the identity matrix. The 
lemma justifies this process and we always assume F is bounded. 

The following theorem is the fundamental result of [23] and gives the required tools to 
address the accuracy of the Riemann-Hilbert numerical methods for orthogonal polynomials 
for arbitrarily large n: 

Theorem 3 Assume 

• C[HI,VL^Y^ exists and the norm ||C[iJ^, ||£(i2(f^j)) < C for all j and n, 

• ||if^||vi/fc,oo(f7j) < C for all j and n, and 

• a{ — )■ as n —7- oo. 

Then for n sufficiently large 
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• The algorithm of Lemma has solutions at each stage, 

• The approximation U^ .^. of U^, the solution of the SIE at stage j in the algorithm of 
Lemma converges in norm, uniformly in n provided mj — )■ oo for all i < j. 

The theorem states that if the contours all have decaying measure then local bounded- 
ness properties on each of the contours can be made global for n large. As we will see below, 
bounding the W'''°° norms of the matrices is often straightforward and the boundedness 
properties of the inverse operator follows from the asymptotic analysis of the RH problem. 

Remark If r„ = fl\ consists of just one scaled contour then the restriction that — )■ 
can be removed due to the fact that z = a\k + is a conformal change of variables for the 
whole problem and this leaves the Cauchy integral operators invariant. 

Remark Similar results hold when the bounds in Theorem |3] are known for a 'nearby' 
RH problem. In this case bounds on the nearby RH problem give slightly weaker convergence 
properties that can still be seen to be uniform in an appropriate sense [23] . 



5.2 The classical Airy parametrix 

In this section we present the deformation and asymptotic solution of the RH problem that 
is performed in the asymptotic analysis of the RH problem for orthogonal polynomials. The 
results from this section can be found in [6]. For brevity of presentation in this section we 
deal with potentials of the form V{x) = x^™. For the asymptotic analysis and deformations 
in the more case of V{x) polynomial see [IHl [HI [HI [9]. 

A sectionally analytic, matrix-valued function $ is constructed explicitly out of the Airy 
function Ai (s) and its derivative such that T<l>^^ — )■ / as n — ?■ oo where T is the solution 
of the original but deformed RH problem. The RH problem for the error E = T<l>^^ has 
smooth solutions and is a near identity RH problem in the sense that the associated singular 
integral operator is expressed in the form I — Kn with ||i^'„||i2 — > as — t- oo. Thus E can 
be computed via a Neumann series for sufficiently large n. 

The deformation proceeds much in the same way as Section 



3.4 



except we replace Pa 

and Pi) with new functions ipa and ipf, that are constructed out of the Airy function. We now 
construct these functions. As an intermediate step, define 



Ai(s) 


Ai (oo^s) 


Ai'(s) 


w^Ai (oj'^s) 


Ai(s) 


Ai (coh) 


Ai'(s) 




Ai(s) 


Ai (coh) 


Ai'(s) 




Ai(s) 


-w^Ai {ojs 


Ai'(s) 


-Ai'ius) 



2-iri 

e 3 . 



1 

-1 1 
1 

1 1 



< arg s < ^ 



27r 
3 



< arg s < IT 



71 < arg s < 



47r 
3 



An 
3 



< args < 27r 



The relations 



Ai (s) + uAi (us) + u^Ai (w^s) 
Ai'(s) + u'^Ai'{us) + uAi'{uj'^s) 
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can be used to show that satisfies the following jump conditions 



^+(s) = ^-(s) < 



See Figure [9] for 7j, i = 1, . . . 4. 





Figure 9: The jump contours for \l/ with jump matrices. We include 6' > in the figure for 
concreteness but its exact value is not needed below. 

Since we only consider V even in this section, the equilibrium measure is supported on a 
symmetric interval [—a, a] for a > 0. Define 



A(z) = -^{z)iz - a)-^'\ \{z) = {z- a)(A(^))2/^ 



(V{z)-i)-g{z). 



It follows from the branching properties of if that A and A are analytic in a neighbourhood 
of a. Furthermore, since A(a) = and A'(a) = (A(a))^/'^ 7^ we use it as a conformal change 
of variables mapping a neighbourhood of z = a into a neighbourhood of the origin. More 
precisely, fix an e > and define Oa = X^^{{\z\ < e}). 
Define 

^aiz) = L(^)^(r^2/=^A(2))e"'^(")"^ 



Liz) 



ipa solves the local RH problem shown in Figure 10(b) The symmetry of V{x) implies that 



^jj-aiz) = a3'ijja{—z)a3 satisfies the jumps shown in Figure 10(a) We are ready to define the 
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full parametrix 



( tpaiz) ZeOa 
<I(Z) = <^ ij.aiz) Z e -Oa 

y N{z) otherwise 




Figure 10: The local parametrices near z = a, —a. As above 6'' > is included for concrete- 
ness but its exact value is not needed, (a) The jump contours for ip_a with jump matrices, 
(b) The jump contours for ipa with jump matrices. 

We will need a result concerning the asymptotics of the Airy function 



2^ \ \s^l^ 

Ai'(s) = --^s3/Vi^'''fl + 0^ ^ 



as s — )■ oo and | args| < vr. These asymptotics, along with the definition of \{z), can be used 
to show 

i;a{z)N~\z) = I + 0{n-^), z e dOa, (6) 
^p.a{z)N--\z) = I + 0{n-^), zedO-a, (7) 

as n — )■ oo uniformly in z provided Oa U 0_a is contained in a sufficiently narrow strip 
containing the real line. See [6] for the details. 

We take the RH problem for T in Figure |6] and label dOa and dO-a- Note that without 
loss of generality we take Oa and 0_a to be open balls around a and —a, respectively. 
Analyticity allows us to deform any open, simply connected set containing a or —a to a ball. 

Since ipa and ip^a solve the RH problem locally in Oa and 0_a, respectively, the function 



E = is analytic in Oa and 0-a- See Figure [llj for the jump contour, fi, and jump 

matrix, J, for the RH problem for E. It is shown in [6] using ^ that the jump matrix for 
this RH problem tends uniformly to the identity matrix as n — )■ oo, again provided that all 
contours are in sufficiently small neighbourhood of the real line. Thus 

\\I-C[J; Q] \\c{LHn)) = 0{n~'), 

and a Neumann series will produce the unique solution u of C[J;Q]u = J — I. T is found 
via the expression 

T{z) = {I + Cnu{z))^{z). 
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Figure 11: The jump contours f2 for the error E. The jump matrix J for E which is taken 
as the piecewise definition as shown. 

5.3 Obtaining the bounds in Theorem [3] 

To apply Theorem |3] one has to first identify the correct scahngs for the contours and second, 
estabhsh bounds on the relevant operator norms and function derivatives. 

5.3.1 The RH problem for E 

In this case, we consider numerically solving the RH problem for E, rather than scaling and 
shifting the contours as we do in practice. This simplifies the proof of uniform approximation 
considerably, at the expense of no longer allowing for degenerate potentials, and requiring 
significantly more knowledge in the construction of the RH problem. 

Take r„ = Q; that is, we do not scale the contour. The near-identity nature of the RH 
problem allows us to avoid any scaling of the problem. Using the asymptotic expansions for 
the derivatives of Airy functions one can show that 



-1^ 



II J — /||p^fc,oo(Q)n/^fe(n) — 0{n 

Furthermore, the fact that ||C [J; Q]^^ \\c{L^(n)) < C follows easily from the Neumann series 
argument already given. Thus we expect the numerical method to uniformly approximate 
solutions of this RH problem for small and arbitrarily large n. 

To demonstrate the convergence properties of the solution for large n we use the following 
procedure. Let Um denote the approximation of u obtained using the numerical method for 
RH problems discussed above with m = |m| collocation points per contour. When we break 
up Q into both its non-self-intersecting components and components that can be represented 
by affine transformations of the unit interval we end up with 14 contours. Thus, we use 
a total of 14m collocation points. We solve the RH problem with m = 10 and then again 
with m = 20. We sample Uio at each collocation point for U20 and measure the maximum 
difference at these collocations points. We define this difference to be the Cauchy error. 



Figure 12 demonstrates that the error decreases as n — > 00. 
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1 10 100 1000 lO'* 10^ 10*^ 



n 

Figure 12: The Cauchy error between Uio and U20 as n — )• 00. This plot indicates that it 
takes fewer collocation points to approximate E as n increases. 



5.3.2 The RH problem for $ 

The RH problem that we use in practice, $, is of a fundamentally simpler form. No additional 
special functions (e.g Airy functions) are needed and yet the contours are located away 
from the stationary points, a and h (we return here to allowing general potentials). All 
deformations are performed by a reordering and analytic continuation of previously defined 
functions. 

Assume we are in the non-degenerate case. Represent 

for 

Unfortunately, we have an issue with the jumps on these scaled contours: as n — )■ 00, 
they approach the unbounded singularities of N(z), violating the conditions of Theorem |3} 
However, we can expand 

2/3^ n^/^ fb-aV^ fl -i\ n'^l^ (h-aY(\ \\ ^. 1, 



Letting 



we observe that 



N{a + zn~'^'^)N-^^ and A^a,nA^(a + ^n'^/a)-! 

are uniformly bounded for z restricted to an annulus around zero as — )■ cxd. 
We thus remove the growth in the jumps by conjugating: let 

outside Oa (i.e., the simply connected region surrounding a) and 

$1 = Na,nQl 
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inside Oa- The jumps on dOa thus become: 
and, on the rest of Q"^, 



SO that Qi has bounded jumps. 

Once Qi and thence $i are calculated, we need to bound the jump of $2; which is 

Similar to before, we can find the two-term expansion of N{h+zn~'^/^) to find iVf, „ to perform 
a second conjugation. The asymptotic convergence of $i to / near h ensures that iVfe.„ and 
$1 asymptotically commute. 

We can now appeal to Theorem |3] to show asymptotic stability. The boundedness of the 
jumps can be verified directly. To bound the inverse operators, we can use the boundedness 
of the parametrix of Section 5.2 when the potential is non-degenerate. In the degenerate 
case, we we would need to use the analysis of the RH problem in [Ij. We omit the details 
here for brevity. 



6 Conclusion 

We presented a numerical method for computing statistics of unitary invariant ensembles, 
based on solving the associated Riemann-Hilbert problem numerically. This required solving 
a nonlinear scalar Riemann-Hilbert problem to calculate the g function associated with the 
equilibrium measure. Scaling the contours appropriately resulted in a numerical method 
that remains accurate for large n, without knowledge of the local parametrices. 

Our hope is that this framework will lead to a better understanding of the relationship 
between the potential V , universality laws and finite n statistics. 

A Computing a Hastings— McLeod Solution of the Painleve 
II transcendent 

Here we focus on the (homogeneous) Painleve II ODE, it is as follows: 

u"{x) = xu{x) + 2u^{x). (8) 

(For brevity we refer to the homogeneous Painleve II simply as Painleve II.) There are many 
important applications of this equation: the Tracy-Widom distribution ^\ from random 
matrix theory is written in terms of the Hastings-McLeod solution [12] and asymptotic 
solutions to the Korteweg-de Vries and modified Korteweg-de Vries equations can be written 
in terms of Ablowitz-Segur solutions [Ij . The aim of this section is to demonstrate that the 
RH formulation can indeed be used effectively to compute solutions to Painleve II, even in 
the asymptotic regime. 

Solutions to differential equations such as ([s]) are typically defined by initial conditions: 
at a point x we are given u{x) and u'{x). In the RH formulation, however, we do not 
specify initial conditions. Rather, the solution is specified by the Stokes' constants] constants 
si, S2, S3 which satisfy the following condition: 
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Assumption 3 



Sl- S2 + S3 + S1S2S3 = 0. 



(9) 



We will treat the Stokes' constants as given, as, in many applications they arise naturally 
whilst initial conditions do not. Given such constants, we denote the associated solution to 
Pby 



Pn{si,S2,S3;z). 



(10) 



Pu and its derivative can be viewed as the special function which map Stokes' constants to 
initial conditions. 

At first glance, computing solutions to (|8| appears trivial: given initial conditions, simply 
use one's favorite time-stepping algorithm, or better yet, input it into an ODE toolbox 
such as Matlab's ode45 or Mathematica's NDSolve. Unfortunately, several difficulties 



immediately become apparent. In Figure 13, we plot several solutions to (|8]) (computed 
using the approach we are advocating): the Hastings-McLeod solution and perturbations 
of the Hastings-McLeod solution. Note that the solution is inherently unstable, and small 
perturbations cause oscillations — which make standard ODE solvers inefficient — and poles 
— which will completely break such ODE solvers (though this issue can be resolved using 
the methodology of [S]). 





(a) 



(b) 



Figure 13: Solutions to Painleve II. (a) Radically different solutions for x < 0. (b) Radically 
different solutions for x > 0. 

Remark There are many other methods for computing the Tracy-Widom distribution 
itself as well as the Hastings-McLeod solution [3l [2], based on the Fredholm determinant 
formulation or solving a boundary value problem. Moreover, accurate data values have been 
tabulated using high precision arithmetic with a Taylor series method |2H However, 
we will see that there is a whole family of solutions to Painleve II which exhibit similar 
sensitivity to initial conditions, and thus a reliable, general numerical method is needed even 
for this case. 

Let $(a;; A) solve the RH problem depicted in Figure 14 let F = Fi U ■ ■ ■ U Fg for 
P^|^gm(«/3-i/6) . g g ^+1^ ^ Y consists of six rays emanating from the origin, as see in 



Figure 14 Then the jump matrix is defined by G{x; A) = G^ix; A) for z G F^, where 



G^{x;\) = G^i\) 



1 „ ^-i8/3\^-2ixX 

1 

1 

„ ^i8/3\^+2ixX 1 



if K, even, 
if K odd. 
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Figure 14: The contour and jump matrix for the Painleve II RH problem. 



This is the RH problem which was solved numerically in [2T]. We can recover the corre- 
sponding solution to Painleve II from $ by [r2j 

Pii{si, S2, S3; x) = 2 lim X^{x; A)i2. 

A— >oo 

As |x| becomes large, the jump matrices G are increasingly oscillatory. We will combat 
this issue by deforming the contour so that these oscillations be exponential decay. To 
simplify this procedure, we first rescale the RH problem. Note that, if we let z = y^\x\X, 
then the jump contour F remains unchanged, and 

<l>+(^) = $+(v^A) = $-(v1^A)G(v^A) = ^-{z)G{z), 

where G{z) = G^iz) on F^ for 



GJz) 



1 s^e-^l^l''"^(^) 
1 

1 



if K even, 



if K odd, 



and 



Then 



= - (Az^ + 2e''''^''z) . 



S2, S3; x) = 2i lim A$(x; A)i2 = 2i\^ lim z)i2. 

A— >oo A— >oo 



A.l Positive x with si = 

We will now deform the RH problem for Painleve II so that numerics are asymptotically 
stable for positive x. We will see that the deformation is extremely simple under the following 
assumption: 
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Assumption 4 S2 = 



We remark that, unlike other deformations, the following deformation can be easily ex- 
tended to achieve asymptotic stability for x in the complex plane such that — | < arga; < |. 

On the undeformed contour, the terms e.^AA^^^s(z) i^g^ome oscillatory as |a;| becomes large. 
However, with the right choice of curve h{t)^ ^±te{h{t)) j-^^^g oscillations; instead, it decays 
exponentially fast as t — )■ oo. But h is precisely the path of steepest descent, which passes 
through the stationary points of 9, i.e., the points where the derivative of 9 vanishes. We 
readily find that 

9'{z) = 2{Az' + l), 

and the stationary points are z = ±z/2. 

We note that, since G2 = I, when we deform Fi and through i/2 they become 
completely disjoint from r4 and Fg, which we then deform through —i/2. We also point out 
that 6*3 ^ = Gi and Gq^ = G4,; thus we can reverse the orientation of and resulting in 



the jump Gi on the curve Ff and G4 on F^^, as seen in Figure 15 




Gi 



G3 



Gi 




Gi 



Ge, 
Gi 



Figure 15: Deforming the RH problem for positive x, with Assumption |4} 
Now recall that 



^ ± 



I 

±-. 
3 



However, we only have F^- emanating from z/2, with jump matrix 



1 



This is exponentially decaying to the identity along F^^; as is G4 along F^^. We will employ the 
approach of Section |5} We first use Lemma |4] to truncate the contours near the stationary 
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point. What remains is to determine what near means. Because 6 behaves hke (9 (z ± i/^Y 
near the stationary points, Assumption [T] imphes that we should choose the shifting of 
/3i = i/2 and (32 = the scahngs ai = a2 = r\x\~^^^ and the canonical domains 

Qi = Q2 = [—1, !]• Here r is chosen so that what is truncated is negligible in the sense of 
Lemma |4} Gq is similar. The complete proof of asymptotic stability of the numerical method 



proceeds in a similar way as in Section 5.3.1 

A. 2 Negative x with si = —S3 = zti and S2 = 

We now develop deformations for the Hastings-McLeod solution for negative x, which cor- 
responds to si = ±i, S2 = and S3 = =Fi [12\- We will realize numerical asymptotic stability 
in the aforementioned sense. 

Assumption 5 si = —S3 = ±i and S2 = 



We begin by deforming the RH problem (Figure 14) to the one shown in Figure 16 The 



horizontal contour extends from —a to a for a > 0. We will determine a below. Define 

-j|x|3/26»(z) 



Co — G()Gi 



Sie 



j|a:|3/26l(z) 



Note that the assumption S2 = simplifies the form of the RH problem substantially, see 



Figure 16(b) We use an approach similar to that of the equilibrium measure to replace 6 






G2 




GeGi _ ^ 


G5 






Go 




(a) 



(b) 



Figure 16: Deforming the RH problem for negative x, with Assumption [5j The black dots 
represent ±a. (a) Initial deformation, (b) Simplification stemming from Assumption [s] 

with a function possessing more desirable properties. Define 

1 



0(z) = e ' ' 2 , CTs 



The branch cut for g{z) is chosen along [—a, a\. If we set a = l/\/2 the branch of g can be 
chosen so that (2;)—^ (z) ~ 0{z~^). Furthermore, 5^+ (2) + (2;) = Q axidlTii{g^{z)—g^{z)) > 
on {—a, a). Define G,- = BI^C-B-l and note that 



G.iz] 



j|^|3/2 g+(^)+g-W 
., |3/2 9+W+9ziW ,-|^|3/2 g-(^)-9-(^) 



Sie 



si e 



Sl 

;|j.|3/2 g-(^)-3-(^) 
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As a; — )■ — oo, Go tends to the matrix 



J 



Sl 



The solution of the RH problem 

^+(2) = ^-(2)J, ze[-a,a], ^(00)=/, 



is given by 



z — a 



z + a 



1/4 



Here (3 has a branch cut on [—a, a] and satisfies I3{z) — )■ 1 as 2 — )■ 00. It is clear that 
-,out^ ^ J uniformly on every closed sub interval of {—a, a). 

We define local parametrices near ±a: 



3 < arg(z — a) < I 



HM 



I if 

G^f ^ if f < a,ig(z — a) < TT 

Gi if — TT < arg(z — a) < 

/ if ^ 



HM 



y-1 

Gi if 



< arg(2; + a) < vr or — tt < arg(2; + a) < 

27r 



2n 
' 3 



Gi if < arg(2: + a) < 3 

^ < arg(z + a) < ^ 



We are ready to define the global parametrix. Given r > define 







if 


z — a 


< r 






if 


z + a 


< r 






if 


z + a 


> r and 2; — a > r 



It follows that \i/HM satisfies the RH problem shown in Figure A.2 

\ 

Figure 17: The jump contours and jump matrices for the RH problem solved by ^E'hm- The 
radius for the two circles is r. 




Let $ be the solution of the RH problem shown in Figure 



18 (a: 



It follows that A = $\I^HM 



solves the RH problem shown in Figure 18(b) The RH problem for A has jump matrices 
that decay to the identity away from ±a. We use Assumption [T] to determine that we should 
We solve the RH problem for A numerically. To compute the solution of 



use r 



\x\ 



Painleve II we use the formula 

Pii{±i,0,^i;x) = 2i lim zA{z)i2. 
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Figure 18: The final deformation of the RH problem for negative x, with Assumption |5j The 
black dots represent ±a. (a) After conjugation by 9. (b) Bounding the contours away from 
the singularities of g and /3 using \1/jjm- 



See Figure 19(a) for a plot of the Hastings-McLeod solution with si = i. To verify our 



computations we may we use the asymptotics 



Pii(i,0, -i;x) ~ 



-X 



'111 



We define 



Pii(z,0, -i]x) + 



X 



-5/2 



to be the relative error which should tend to a constant for x large and negative. We 



demonstrate this in Figure 19(b) 



Remark Since (3 has unbounded singularities we expect that a similar issue as in Sec- 



tion |5.3.2| will arise. We do not go though the details of this but this approach produces 
accurate numerics for all x on the real line. 
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